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We apply the Lagrange Multiplier method to study the uncertainties of physical 
predictions due to the uncertainties of parton distribution functions (PDFs), using the cross 
section aw for W production at a hadron collider as an archetypal example. An effective \ 2 
function based on the CTEQ global QCD analysis is used to generate a series of PDFs, each of 
which represents the best fit to the global data for some specified value of aw- By analyzing 
the likelihood of these "alterative hypotheses", using available information on errors from 
the individual experiments, we estimate that the fractional uncertainty of aw due to current 
experimental input to the PDF analysis is approximately ±4% at the Tevatron, and ±8-10% 
at the LHC We give sets of PDFs corresponding to these up and down variations of aw- We 
also present similar results on Z production at the colliders. Our method can be applied to 
any combination of physical variables in precision QCD phenomenology, and it can be used 
to generate benchmarks for testing the accuracy of approximate methods based on the error 
matrix. 
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1 Introduction 



All calculations of high energy processes with initial hadrons, whether within the Standard 
Model (SM) or exploring New Physics, require parton distribution functions (PDFs) as an 
essential input. The reliability of these calculations, which underpins both future theoretical 
and experimental progress, depends on understanding the uncertainties of the PDFs. The 
assessment of PDF uncertainties has, therefore, become an important challenge to high 
energy physics in recent years. 

The PDFs are derived from global analyses of experimental data from a wide range of 
hard processes in the framework of perturbative quantum chromodynamics (PQCD). Quan- 
tifying the uncertainties in a global QCD analysis is far from being a straightforward exercise 
in statistics. There are non-gaussian sources of uncertainty from perturbation theory (e.g., 
higher order and power law corrections), from choices of parametrization of the nonpertur- 
bative input (i.e., initial parton distributions at a low energy scale), from uncertain nuclear 
corrections to experiments performed on nuclear targets, and from normal experimental sta- 
tistical and systematic errors. These sources of error need to be studied individually, and 
eventually combined in a systematic way. 

We shall be concerned in this paper with uncertainties of PQCD predictions due to 
uncertainties of PDFs arising from experimental measurement errors. This problem is con- 
siderably more complicated than it appears on the surface. The reason is that in a global 
analysis, the large number of data points (~ 1300 in our case) do not come from a uniform set 
of measurements, but consist of a collection of measurements from many experiments (~ 15) 
on a variety of physical processes (~ 5 — 6) with diverse characteristics, precision, and error 
determination. The difficulty is compounded by a large number of fitting parameters (~ 16) 
which are not uniquely specified by the theory. Several approaches to this problem have been 
proposed, with rather different emphases on the rigor of the statistical method, scope of ex- 
perimental input, and attention to various practical complications jl], 0, ^, ffl, [5|, ^ |7J . Our 
group has initiated one of these efforts, with the emphasis on utilizing the full constraints of 
the global data 0. This work has motivated a closer examination of the standard techniques 
of error analysis, and necessary improvements and extensions to these techniques, as applied 
to a complex real world problem such as global QCD analysis of PDFs |§. 

In this paper we present a detailed analysis of uncertainties of physical observables 
due to parton distribution functions, using the Lagrange Multiplier method proposed in 
|7], §]. This method explores the entire multi-dimensional parton parameter space, using 
an effective x 2 function that conveniently combines the global experimental, theoretical, 
and phenomenological inputs to give a quantitative measure of the goodness-of-fit for a 
given set of PDF parameters. (Cf. Sec. 0.) The method probes directly the variation of 
the effective \ 2 along a specific direction in the PDF parameter space — that of maximum 
variation of a specified physical variable. The result is a robust set of optimized sample 
PDFs (or "alternative hypotheses") from which the uncertainty of the physical variable 
can be assessed quantitatively without the approximations inherent in the traditional error 
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matrix approach. For concreteness, we consider the cross section aw of W boson production 
at the Tevatron as the archetypal example. {Cf. Sec. |[) 

The definition of the effective \ 2 function, and the inputs that go into it, do not permit a 
direct statistical interpretation of its numerical value. To obtain meaningful confidence levels 
for the optimized sample PDF sets, it is necessary to conduct a series of likelihood analyses of 
these sample PDFs, using all available information on errors for the individual experiments. 
The results from these analyses serve as the basis to assign an overall uncertainty range on 
the physical variable, and a corresponding tolerance measure for the effective \ 2 function 
used in the analysis, that are consistent with the experiments used in the current global 
QCD analysis. {Cf. Sec. |.) 

This method can be applied to any physical variable, or combination of physical vari- 
ables, in precision QCD phenomenology. In Sec. |5] we present results on W production at 
the LHC, and Z production at the Tevatron and the LHC We compare the uncertainties 
obtained in all cases, and comment on previous estimates in the context of these results. In 
Sec. H we present parton distribution sets that are optimized to give high/low values of the 
W and Z cross sections, while remaining consistent with current experiments according to 
our analysis. 

The Lagrange Multiplier method provides a useful tool to test the reliability of the more 
traditional method of error propagation via the error matrix |1|, |], |9|, which relies on the 
quadratic expansion of the x 2 function around its minimum. In a companion paper we 
perform an in-depth analysis of the uncertainties of the PDFs in the error matrix approach, 
using the much improved numerical method for calculating the Hessian that was developed 
in || . There we demonstrate how the more specialized Lagrange Multiplier method can set 
useful benchmarks for the general purpose error matrix approach. 



2 The Global QCD Analysis 



We adopt the same experimental and theoretical input as the CTEQ5 analysis |ll[]: 15 data 



sets from 11 experiments on neutral-current and charged- current deep inelastic scattering 
(DIS), lepton-pair production (Drell-Yan), lepton asymmetry in H^-production, and high 
inclusive jet production processes are used. {Cf. Table [I] in Sec. £|.) The total number of data 
points is iV = 1295. We denote the experimental data values by {D} = {Dj; 1 = 1,..., N}. 
The theory input is next-leading-order (NLO) PQCD, and the theory value for the data 
point / will be denoted by Tj. The theory depends on a set of parameters {a} = {af, i = 
1, . . . ,d}. These parameters characterize the nonperturbative QCD input to the analysis; 
they determine the initial PDFs {f{x, Q ; {a})} defined at a low energy scale Q , below the 
energy scale of the data, which we choose to be Qq = 1 GeV. When we need to emphasize 
that the theoretical values depend on the PDF parameters we write Tj{a) to indicate the 
dependence on {a}. 

The parametrization of {f{x, Qq)} is somewhat arbitrary, motivated by physics, nu- 
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merical considerations, and economy. Another parametrization might be employed, and 
differences among the possible parametrizations are in principle a source of theoretical un- 
certainty in their own right. For most of this study we focus on a single parametrization, 
but we comment on the effect of changing the parametrization at the end of Sec. f|. The 
number d of the parameters {a} is chosen to be commensurate with current experimental 
constraints. For this study we use d = 16. The detailed forms adopted for the initial functions 
{f(x, Qo] {«})} are not of particular concern in this study, since we shall be emphasizing re- 
sults obtained by ranging over the full parameter space. f\ The explicit formulas are given in 
Appendix y (where relevant PDFs from the results of our study are presented). The T/({a}) 
are calculated as convolution integrals of the relevant NLO QCD matrix elements and the 
universal parton distributions {f(x,Q; {a})} for all Q. The latter are obtained from the 
initial functions {f(x, Qo] {a})} by NLO QCD evolution. 

The global analysis consists of a systematic way to determine the best values for the 
{a}, and the associated uncertainties, by fitting {T(a)} to {D}. Because of the wide range 
of experimental and theoretical sources of uncertainty mentioned in the Introduction, there 
are a variety of strategies to deal with the complex issues involved ||, 0, £|, [7[]. In the 
next two sections, the primary tool we employ is conventional x 2 analysis. The important 
task is to define an effective x 2 function, called Xgiobai( a )> that conveniently combines the 
theoretical and global experimental inputs, as well as relevant physics considerations based 
on prior knowledge, to give an overall measure of the goodness-of-fit for a given set of PDF 
parameters. 

Experience in global analysis of PDFs during the past two decades has demonstrated 
that the PDFs obtained by the minimization of such a suitably chosen Xgiobai provide very 
useful up-to-date hadron structure functions which, although not unique, are representative 
of good fits between theory and experiments. Now we must quantify the uncertainties of 
the PDFs and their predictions; i.e., we must expand the scope of the work from merely 
identifying typical solutions to systematically mapping the PDF parameter space in the 
neighborhood around the minimum of x 2 - 

The simplest possible choice for the x 2 function would be 

*■(«) = E [J '~? (a)1 ' « 

7=1 a I 

where 07 is the error associated with data point /. Through Tj(a), x 2 ( a ) is a function of the 
theory parameters {a}. Minimization of x 2 { a ) would identify parameter values for which 
the theory fits the data. However, the simple form (0) is appropriate only for the ideal case 
of a uniform data set with uncorrelated errors. For data used in the global analysis, most 
experiments combine various systematic errors into one effective error for each data point, 
along with the statistical error. Then, in addition, the fully correlated normalization error 

a In other words, for this paper, the PDF parameters {a} play mostly the role of "internal variables" . In 
contrast, they occupy the center stage in the companion paper JTc| ] . 
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of the experiment is usually specified separately. For this reason, it is natural to adopt the 
following definition for the effective x 2 ( as done in previous CTEQ analyses): 

Xgiobai( a ) = ^w„Xn( a ) (n labels the different experiments) (2) 

n 

xM = [—) +E{ ^ ) (3) 

For the n th experiment, D n j, a®i, and T n i(a) denote the data value, measurement uncer- 
tainty (statistical and systematic combined), and theoretical value (dependent on {a}) for 
the I th data point; is the experimental normalization uncertainty; M n is an overall nor- 
malization factor (with default value 1) for the data of experiment n. The factor w n is a 
possible weighting factor (with default value 1) which may be necessary to take into account 
prior knowledge based on physics considerations or other information. The a priori choices 
represented by the w n values are present, explicitly or implicitly, in any data analysis. For 
instance, data inclusion or omission (choices which vary for different global analysis efforts) 
represent extreme cases, assigning either 100% or 0% weight to each available experimental 
data set. Similarly, choices of various elements of the analysis procedure itself represent 
subjective input. Subjectivity of this kind also enters into the analysis of systematic errors 
in experiments. 

The function Xgiobai( a ) allows the inclusion of all experimental constraints in a uniform 
manner while allowing flexibility for incorporating other relevant physics input. We will 
make use of this function to explore the neighborhood of the best fit, and to generate 
sample PDFs pertinent to the uncertainty of the prediction of a specific physical variable 
of interest. However, the numerical value of this effective \ 2 function should not be given 
an a priori statistical interpretation, because correlations between measurement errors, and 
correlated theoretical errors, are not included in its definition. In particular, the likelihood of 
a candidate PDF set {a} cannot be determined by the value of the increase AXg lobal (a) above 
the minimum.0 Instead, the evaluation of likelihoods and estimation of global uncertainty 
will be carried out in a separate step in Sec.£|, after sets of optimal sample PDFs for the 
physical variable of interest have been obtained. 



3 The Lagrange Multiplier Method 

The Lagrange Multiplier method is an extension of the x 2 minimization procedure, that 
relates the range of variation of a physical observable X dependent upon the PDFs, to the 
variation of the function Xgiobai( a ) that is used to judge the goodness of fit of the PDFs to 
the experimental data and PQCD. 

b The often quoted theorem of Gaussian error analysis, that an increase of x 2 by 1 unit in a constrained 
fit to data corresponds to 1 standard deviation of the constrained variable, is true only in the absence of 
correlations. When existing correlations are left out, the relevant size of A% 2 can be much larger than 1. 
Appendix H discusses this point in some detail. 
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3.1 The Method 



The method has been introduced in |7], || . The starting point is to perform a global analysis as 
described in Sec. ||, by minimizing the function Xgiobai( a ) defined by Eq. (§), thus generating a 
set of PDFs that represents the best estimate consistent with current experiment and theory. 
We call this set the "standard set"f], denoted S . The parameter values that characterize 
this set will be denoted by {a^} = {af^; i — 1, . . . , d}; and the absolute minimum of Xgiobai 
will be denoted by Xo- Now, let X be a particular physical quantity of interest. It depends 
on the PDFs, X = X(a), and the best estimate (or prediction) of X is Xo = X^a^). We 
will assess the uncertainty of this predicted value by a two-step analysis. First, we use the 
Lagrange Multiplier method to determine how the minimum of Xgiobai( a ) increases, i.e., how 
the quality of the fit to the global data set decreases, as X deviates from the best estimate 
X . Second, in Section 4, we analyze the appropriate tolerance of Xgiobai- 

As explained in [0, |j, the first step is taken by introducing a Lagrange multiplier 
variable A, and minimizing the function 

*(A,a) = XgiobaiW + AX(a) (4) 

with respect to the original d parameters {a} for fixed values of A. In practice we minimize 
\I/(A, a) for many values of the Lagrange multiplier: Ai, A2, . . . , Am- For each specific value 
Aq,, the minimum of ty(\ a ,a) yields a set of parameters {a min (A a )}, for which we evaluate 
the observable X and the related Xgiobar We use the shorthand (X a , Xgiobaia) f° r this pair. 
Xgiobai.a represents the lowest achievable Xgiobai' f° r the global data, for which X has the 
value X a , taking into account all possible PDFs in the full d- dimensional parameter space 
of points {a}. In other words, the result {a m i n (A a )} is a constrained fit — with X constrained 
to be X a . We can equivalently say that X a is an extremum of X if Xgiobai i s constrained to 
he Xgiobai.a- We denote the resulting set of PDFs by S a . 

We repeat the calculation for many values of A, following the chain 

\ a — > min [®(\ a , a)] — > a min (A Q ) — > X a and Xgiobai,a 

for a — 1, 2, 3, ... , M. The result is a parametric relationship between X and Xgiobai> through 
A. We call this function XgiobaiW; so Xgiobai( x «) = Xgiobai,a is th e minimum of XgiobaiO) 
when X is constrained to be X a . The absolute minimum of XgiobaD which we denote xh is 
the minimum of \1/(A = 0, a), occurring at {a} = {a < -°- ) }. Thus the procedure generates a set 
of optimized sample PDFs along the curve of maximum variation of the physical variable X 
in the <i-dimensional PDF parameter space (with d = 16 in our case). These PDF sets {S a } 
are exactly what is needed to assess the range of variation of X allowed by the data. In 
other words, the Lagrange Multiplier method provides optimal PDFs tailored to the physics 
problem at hand, in contrast to an alternative method |§ that generates a large sample 
of PDFs by the Monte Carlo method. The underlying ideas of these two complementary 
approaches are illustrated in the plot on the left side of Fig. [l]. 
c This standard set is very similar to the published CTEQ5M1 set p[~T[ ] . 
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Figure 1: Left: The LM method provides sample points along a single curve Lx in the 
multi-dimensional PDF parameter space, relevant ro the observable X. Right: For a given 
tolerance Axg lobal , the uncertainty in the calculated value of X is ±AX. The solid points 
correspond to the sample points on the curve L x in the left plot. 

Xgiobai(^) * s the lowest achievable value of Xgiobai( a ) f° r the value X of the observable, 
where Xgiobai( a ) represents our measure of the goodness-of-fit to the global data. Therefore 
the allowed range of X, say from Xq — AX to Xo + AX, corresponding to a chosen tolerance 
of the goodness of fit AXg lobal = Xgiobai — Xo> can be determined by examining a graph of 
Xgiobai versus X, as illustrated in the plot on the right side of Fig. [TJ. This method for 
calculating AX may be more robust and reliable than the traditional error propagation 
because it does not approximate X(a) and Xgiobai( a ) by linear and quadratic dependence on 
{a}, respectively, around the minimum. 

Although the parameters {a} do not appear explicitly in this analysis, the results do 
depend, in principle, on the choice of parameter space (including the dimension, d) in which 
the minimization takes place. In practice, if the degrees of freedom represented by the 
parametrization are chosen to match the constraining power of the global data sets used, 
which must be true for a sensible global analysis, the results are quite stable with respect to 
changes in the parametrization choices. The sensitivity to these choices is tested, as part of 
the continuing effort to improve the global analysis. 

The discussion so far has left open this question: What is the appropriate tolerance 
AXg lobal to define the "error" of the prediction X ? This question will be addressed in Sec. [|[ 

Our method can obviously be generalized to study the uncertainties of a collection 
of physical observables (X 1; X 2 , . . . ,X S ) by introducing a separate Lagrange multiplier for 
each observable. Although the principle stays the same, the amount of computational work 
increases dramatically with each additional observable. 
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Figure 2: Calculated cross section for W ± boson production (multiplied by the branching 
ratio for W~ — > eV) at the Tevatron, for various current and historical PDFs. The two plots 
are from Refs. and 0] respectively. 



3.2 A case study: the W cross section 

In this subsection we examine the cross section aw for inclusive W ± production at the 
Tevatron (pp collisions at ^/s = 1.8 TeV) to illustrate the method and to lay the ground 
work for the quantitative study of uncertainties to be given in Sec. 4. Other examples will be 
described in Sec. ||. Preliminary results of this section have been reported previously |7], ||. 

Until recently the only method for assessing the uncertainty of aw due to PDFs has 
been to compare the calculated values obtained from a number of different PDFs, as illus- 
trated in Fig. |2|, in which the plots are taken from existing literature.^ The PDFs used in 



these comparisons are either the "best fits" from different global analysis groups [JTTJ, [12 
(hence are not pertinent to uncertainty studies) or are chosen by some simple intuitive cri- 
teria [ 1~3| . The meaning and reliability of the resulting range of aw are not at all clear. 



Furthermore, these results do not provide any leads on how the uncertainties can be im- 
proved in the future. The Lagrange Multiplier technique provides a systematic method to 
address and remedy both of these problems. 

Let the physical quantity X of the last subsection be the cross section aw for W ± 
production at the Tevatron. Applying the Lagrange method, we obtain the constrained 
minimum of Xgiobai as a function of aw, shown as solid points in Fig. |3|. The best esti- 
mate value, i.e., the prediction for the standard set So, is awo = 21.75 nb. The curve is a 
polynomial fit to the points to provide a smooth representation of the continuous function 
XgiobaiPO- We see that all the sample PDF sets obtained by this method lie on a smooth 
quasi-parabolic curve with the best-fit value at the minimum. 

As discussed earlier (in Fig. |l|) points on the curve represent our sample of optimal 

d These plots show the product of aw times a leptonic branching ratio, which is what is measured exper- 
imentally. The branching ratio B has some experimental error. For studying the uncertainties of aw, we 
will focus on aw itself in the rest of the paper. 
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W production at the Tevatron 
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Figure 3: Minimum Xgiobai versus aw, the inclusive W ± production cross section at the 
Tevatron {pp collisions at s/s = 1.8 TeV) in nb. The points were obtained by the Lagrange 
Multiplier method. The curve is a polynomial fit to the points. 

PDFs relevant to the determination of the uncertainty of aw To quantify this uncertainty, 
we need to reach beyond the effective Xgiobai function, and establish the confidence levels for 
these "alternative hypotheses" with respect to the experimental data sets used in the global 
analysis. 




4 Quantifying the Uncertainty 

Consider a series of sample PDF sets along the curve Xgiobai(^0 °f Fig- H denoted by {S a ; a = 
0, 1, . . . , M} where So is the standard set. These represent "alternative hypotheses" for the 
true PDFs, and we wish to evaluate the likelihoods associated with these alternatives. To 
do so, we go back to the individual experiments and, in each case, perform as detailed a 
statistical analysis as is permitted with available information from that experiment. After we 
have obtained meaningful estimates of the "errors" of these candidate PDFs with respect to 
the individual experiments, we shall try to combine this information into a global uncertainty 
measure in the form of AX and AXg lobal . 

The experimental data sets included in our global analysis are listed in Table |l|. For 
some of these experiments, information on correlated systematic errors is available (albeit 
usually in unpublished form). For these, statistical inference should be drawn from a more 
accurate x\ function than the simple formula Eq. (j^) used for the global fit. In particular, 
if a n i is the uncorrelated error and {Pki] k = 1,2,..., K} are the coefficients of K distinct 
correlated errors associated with the data point /, then an appropriate formula for the Xn 
function is 

a=E ( % r ^ -x;x>fr*-') M ,*» <«) 

I nI k=l k'=l 
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Experiment Process 



Label 



# Data pts Reference 



BCDMS 
BCDMS 
HI 

ZEUS 

NMC 

NMC 

NMC 

CCFR 

CCFR 

E605 

NA51 

E866 

CDF 

DO 

CDF 



DIS HP 
DIS fid 
DIS ep 
DIS ep 
DIS up 
DIS HV/^ n 
DIS hp/ ^ n 
DIS up 
DIS i/p 
D-Y pp 
D- Y pc? / pp 
D-Y pd/pp 

^ T lep—asym. 

pp — > jetX 
pp — > jet X 



BCDMSp 
BCDMSd 
HI 

ZEUS 

NMCp 

NMCr 

NMCrx 

CCFR2 

CCFR3 

E605 

NA51 

E866 

CDFw 

DOjet 

CDFjet 



168 
156 
172 
186 
104 
123 

13 
87 
87 
119 

1 

11 
11 

24 
33 



14 



14 



El 







171 









20 



21 





E3 



24 



Table 1: List of data sets used in the global analysis, 
where is a vector, and A^y a matrix, in K dimensions: 

Bk — y^hljDnl — Tnlj/Pni ] A k k' — £fcfc' +/ J PkI^k'l/o'ni- (6) 
I I 

(The sum over / here includes only the data from experiment n.) Traditionally, Xn is written 
in other ways, e.g., in terms of the inverse of the (N x N) variance matrix. For experiments 
with many data points, the inversion of such large matrices may lead to numerical insta- 
bilities, in addition to being time-consuming. Our formula has a significant advantage 
in that all the systematic errors are first combined ("analytically") in the definitions of Bk 
and Akk'- Equation ([5]) requires only the inverse of the much smaller (K x K) matrix Akw- 
(K is the number of distinct systematic errors.) The derivation of these formulas is given 
in Appendix 0. Equation (H) reduces to the minimum of xi m Eq. (|3]) with respect to Af n 
if the only correlated error is the overall normalization error for the entire data set; in that 
case /3j = -a%D n i- 

By using Eq. ([5|), or Eq. (|^) for cases where the correlations of systematic errors are 
unavailable, we obtain the best estimate on the range of uncertainty permitted by available 
information on each individual experiment. We should note that the experimental data sets 
are continuously evolving. Some data sets in Table 1 will soon be updated (Zeus, HI) or 
replaced (CCFR).0 In addition, most information on correlated systematic errors is either 

e Cf. Talks presented by these collaborations at DIS2000 Workshop on Deep Inelastic Scattering and 
Related Topics, Liverpool, England, April 2000. 
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unpublished or preliminary. The results presented in the following analysis should therefore 
be considered more as a demonstration of principle — as the first application of our proposed 
method — rather than the final word on the PDF uncertainty of the W cross section. 

4.1 Uncertainty with respect to individual experiments 

As an example, we begin by comparing the {S a } series for aw at the Tevatron to the HI 
data set fIS| . Results on correlated systematic errors are available for this data set ,[| and 
are incorporated in the calculation using Eq. (|5|). The number of data points in this set 
is N H1 = 172. The calculated values of Xm/Nm are plotted against aw in Fig. [|. The 
curve is a smooth interpolation of the points. The value of Xhi/Nhi for the standard set So 
(indicated by a short arrow on the plot) is 0.975; and it is 0.970 at the minimum of the curve. 
These values are quite normal for data with accurately determined measurement errors. We 
can therefore apply standard statistics to calculate the 90% confidence level on x 2 /N for 
iV = 172. The result is shown as the dashed horizontal line in Fig. |j. 



1.1 
%1.06 
S 1.02 
0.98 

21 21.5 22 22.5 
cr w (nb) Tevatron 

Figure 4: x 2 /N °f the HI data, including error correlations, for sample PDFs obtained by 
the Lagrange Multiplier method for constrained values of aw at the Tevatron. The arrow 
indicates the global minimum. 

We have similarly calculated Xn/N n including information on the correlations of sys- 
tematic errors for the BCDMSp data set. The results are similar to the HI results, except 
that the absolute values are all larger than 1.12, a large value for iV = 168 data points. This 
is a familiar problem in data analysis, and it is encountered in several other data sets in this 
global analysis (cf. below). The x^/iV n calculation including correlations of the errors is also 
done for the DO and CDF jet cross sections. | For those experiments that have only provided 
(effective) uncorrelated errors, we must rely on Eq. (|3|) for our error calculation, since that 
represents the best information available. 

f These systematic errors are unpublished results, but are made available to the public on the HI Web 
page. For convenience, we have approximated each of the pair of 4 non-symmetrical errors by a single 
symmetric error. The size of the resulting error on aw inferred from this evaluation is not affected by that 
approximation. 

g The measurement errors of the jet cross sections are dominated by systematic errors, so the error 
correlation matrices are used for Xn °f these experiments even in Xgiobai- 
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In order to obtain usable likelihood estimates from all the data sets, one must address 
the problem mentioned in the previous paragraph: Even in a "best fit", the values of x 2 
per data point, x^/N n , for individual experiments vary considerably among the established 
experiments (labeled by n). Specifically, Xn/N n ranges from 1.5 — 1.7 (for ZEUS and CDFjet) 
on the high end to 0.5 — 0.7 (for some D-Y experiments) on the low end in all good fits. 
Considering the fact that some of these data sets contain close to 200 points, the range of 
variation is huge from the viewpoint of normal statistics: Experiments with x 2 n fN n deviating 
from 1.0 by a few times ^J2/N n in either direction would have to be ruled out as extremely 
unlikely f2§. 

The reasons for Xn/N n to deviate from 1.0 in real experiments are complex, and vary 
among experiments. They are, almost by definition, not understood, since otherwise the 
errors would have been corrected and the resulting x 2 would become consistent with the 
expected statistical value. Under these circumstances, a commonly adopted pragmatic ap- 
proach is to focus on the relative x 2 values with respect to some appropriate reference x 2 -Q 
Accordingly, in the context of performing global QCD analysis, we adopt the following pro- 
cedure. For each experiment (labeled by n): 

(i) Let Xno denote the value of xl, f° r t ne standard set So. We assume So is a viable 
reference set. Because Xno ma y be far from a likely value for random errors, we rescale the 
values of x\a. (f° r « = 0, 1, 2, . . . , M) by a factor C n0 , calling the result xi,a 

Xn,a > Xn,a = ^nOX n ,a- CO 

The constant C n0 is chosen such that, for the standard set, Xno assumes the most prob- 
able value for a chi-squared variable: xio = £so = the 50 th percentile of the chi-squared 
distribution P(x 2 , N n ) with N n degrees of freedom, defined by 

I 50 P(x 2 ,N n )d X 2 = 0.50. (8) 
Jo 

(If N n is large then £ 50 « N n .) The rescaling constant C n0 is thus ^5o/Xno- For random 
errors the probability that x 2 < £50 ( or > £50) is 50%. For those experiments whose xl,o 
deviates significantly from £ 50 , this rescaling procedure is meant to provide a simple (but 
crude) way to correct for the unknown correlations or unusual fluctuations. 

(ii) We then examine the values of Xn, a 101 t ne alternative sets S a with a = 1,2, ... , M, 
using Xna ~xio to compute the statistical likelihood of the alternative hypothesis S a with 
respect to the data set n, based on the chi-squared distribution with N n data points. 

h The alternative is to take the absolute values of Xn seriously, and hence only work with alternative 
hypotheses and experiments that are both self-consistent (i.e., have \Xn/N n — 1 1 y/2/N n ) and mutually 
compatible in the strict statistical sense (i.e., have overlapping likelihood functions). Since few of the 
precision DIS experiments are compatible in this sense, one must then abandon global analysis, and work 
instead with several distinct (and mutually exclusive) analyses based on different experiments. 
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This procedure does not affect the results presented earlier for the HI experiment, since 
Xnfl/N n is already very close to 1 for that experiment. 

Before presenting the results of the likelihood calculation, it is interesting to examine, 
in Fig. [|, the differences ^-Xna = Xn, a ~Xn,o (before rescaling) versus aw for the 15 data sets. 
(N.B. The vertical scales of the various plots are not the same, due to the large variations in 
the value of a for different experiments.) The ordering of the experiments in Fig. |5] is the 
same as in Table [I], with experiments ordered by process (DIS, DY, W and jet production). 
It is clear from these graphs that the DIS experiments place the strongest constraints on 
aw, because they have the largest Ax\ for the same Aaw- This is to be expected since 
quark- ant iquark annihilation makes the dominant contribution to aw- We also observe that 
most experiments place some constraint on aw on both sides, but a few bound it on one side 
only. Globally, as shown in Fig. |3], the combined constraints give rise to a classic parabolic 
behavior for Xgiobai^w)- 

X 2 -Xo 2 vscr w (Tevatron) 
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Figure 5: The abcissa is aw in nb, at the Tevatron. The ordinate is Xn ~ Xno- The number 
in parentheses is the number of data points. The horizontal lines are explained in the text. 

To estimate the statistical significance of the individual Xn increases, we assume that 
the rescaled variable Xn obeys a chi-squared distribution P(x 2 , N n ) for N n data points. 
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Thereby, we estimate the value of xl, that corresponds to the 90% confidence level (CL) 
uncertainty for aw (with respect to experiment n) from the formula xl, — £90, where £90 is 
the 90 th percentile defined by 

/ P(x\N n )d X 2 = 0.90. (9) 
Jo 

For example, Fig. || shows the chi-squared distribution P(x 2 , N n ) for N n = 172, the number 
of data points in the HI data set. The 50 th and 90 th percentiles are indicated. We choose a 
conservative 90% CL because there are other theoretical and phenomenological uncertainties 
not taken into account by this analysis. 

0.03 
0.02 
0.01 

130 160 190 220 A 

Figure 6: The chi-square distribution P(x 2 , N n ) for N n = 172 data points. The dashed lines 
indicate the 50 th and 90 th percentiles. 




To summarize our procedure, an alternative PDF set S a lies within the 90% CL for 
experiment n if it has < £90; that is, if 

^ < p-. (10) 



X n ,0 6 



We judge the likelihood of S a from the ratio of x 2 na t° the reference value x\ o> rather than 
from the absolute magnitude. The horizontal lines in Fig. |5] correspond to the values of A-Xn 
obtained in this way. Finally, from the intercepts of the line with the interpolating curve in 
each plot in Fig. [5], we obtain an estimated uncertainty range of aw from each individual 
experiment. The results are presented collectively in Fig. [7], where, for each experiment, the 
point (•) is the value of aw for which x\ i s minimum, and the error bar extends across the 
90% CL based on that data set. 

The uncertainty ranges shown in Fig. [7] with respect to individual experiments represent 
the most definitive results of our study, in the sense that the input and the assumptions can 
be stated clearly and the analysis is quantitative within the stated framework. It is natural 
to proceed further and estimate a global measure of Aaw and the corresponding Ax^ obal . 
This last step is, however, less well-defined and requires some subjective judgement. 
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Figure 7: Ranges of aw within the 90% CL for the individual experiments. The ordinate is 
aw for the Tevatron process pp — > W ± X. The solid line is the best estimate according to 
the standard PDF set S . The dashed lines are the bounds described in Sec. O. 



4.2 The Global Uncertainty 

It should be emphasized that the ranges shown by the error bars in Fig. |7] are not errors 
determined independently by each experiment; rather they represent the ranges allowed by 
the experiments for alternative global fits {S Q }. For this reason, and others related to the 
rescaling of x 2 mentioned earlier as well as approximations inherent in many of the original 
published errors,[] it is not obvious how to combine these errors. We refer to the ranges 
in Fig. by the generic term local (i.e., single-experiment) uncertainties. On a qualitative 
level, Fig. [?| exhibits the same features seen earlier in Fig. |5|: (i) the quark dominated DIS 
experiments give the smallest error bars; and (ii) a few experiments only set bounds on 
one side, while the rest limit the range in both directions. In addition, Fig. |7| gives us an 
overall view which clearly shows that aw is well constrained in the global analysis, and the 
experimental bounds are consistent with each other. 

The important question is how to provide a sensible measure of the overall uncertainty 
in view of the complexity of the problem already described. The situation here is not unlike 
the problem of assigning an overall systematic error to an experimental measurement. Figure 
[7| shows a set of 90% CL ranges for aw from different sources, but these ranges are highly 
correlated, because the alternative hypotheses being tested come from global fits. The final 
uncertainty must be a reasonable intersection of these ranges. 

We will state an algorithm for obtaining the final uncertainty measure of aw based on 
Fig.0. The same algorithm can be applied in the future for predictions of other observables. 

'For instance, the single uncorrelated systematic error associated with each data point, which is the 
only systematic error given for most experimental data sets, is clearly only an "effective uncorrelated error" 
which qualitatively represents the effects of the many sources of systematic error, some of which are really 
correlated. 
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It has two parts: (1) Determine the central value using all the experiments; that is the 
solid line in Fig. (2) Then take the intersection of the error ranges as the combined 
uncertainty. But in calculating the intersection, experiments below the mean are used only 
for setting the lower bound, and experiments above the mean are used only for setting the 
upper bound. With this algorithm, experiments that permit a large range of aw, i-e-, that 
depend on aspects of the PDFs that are not sensitive to the value of aw, will not affect 
the final uncertainty measure (as they should not). According to this algorithm, the result 
for the uncertainty of aw is 20.9nb<o"vK < 22.6 nb. These bounds are approximately ±4% 
deviations from the prediction (21.75 nb) and so we quote a ±4% uncertainty in aw due to 
PDFs. 

Now we may determine the increase in Xgiobai that corresponds to our estimated un- 
certainty Aaw in the aw prediction. Referring to Fig. ||, a deviation of aw by ±4% from 
the minimum corresponds to an increase Axg lobal « 180. That is, Axg lobal in Fig. |l]is 180. In 
other words, along the direction of maximum variation of aw a PDF set with AXg lobal > 180 
is found to violate some experimental constraints by this analysis. 

4.3 Comments 

We should point out that the above uncertainty estimate, Aaw /^w ~ 4%, represents only 
a lower bound on the true uncertainty, since many other sources of error have not yet been 
included in the analysis: theoretical ones such as QCD higher order and resummation effects, 
power-law corrections, and nuclear corrections. These need to be taken into consideration in 
a full investigation of the uncertainties, but that goes beyond the scope of this paper.Q We 
shall add only two remarks which are more directly related to our analysis. 

The first concerns a technical detail. In the results reported so far, we have fixed the 
normalization factors {J\f n } in the definition of Xgiobai (Eq- dD) a t their values determined 
in the standard fit Sq. If we let these factors float when we perform the Lagrange Multiplier 
analysis, Aaw will increase noticeably compared to Fig. |3] for the same A%g lobal . However, 
upon closer examination, this behavior can be easily understood and it does not imply a real 
increase in the uncertainty of aw- The key observation is that the additional increase (or 
decrease) in aw is entirely due to a uniform increase (or decrease) of {Af n } for all the DIS 
experiments. There is a simple reason for this: The values of the q and q distributions in 
the relevant x range (which determine the value of aw) are approximately proportional to 
{Nti}dis- Although every experiment does have a normalization uncertainty, the probability 
that the normalization factors of all the independent DIS experiments would shift in the 
same direction by the same amount is certainly unlikely. Hence we avoid this artificial effect 
by fixing {Af n } at their "best values" for our study. Allowing the factors {J\f n } to vary 
randomly (within the published experimental normalization uncertainties) would not change 
our estimated value of Aaw significantly. 

^Because there are these additional sources of uncertainty, we have used 90% CL's, rather than 68% CL's, 
to calculate the error. 
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The second remark concerns the choice of parametrization. We have mentioned that 
even the robust Lagrange Multiplier method depends in principle on the choice of the parton 
parameter space, i.e., on the choice of the functional forms used for the nonperturbative 
PDFs at the low momentum scale Qq. To check how our answers depend on the choice of 
parametrization in practice, we have done many similar calculations, using different numbers 
of free parameters within the same functional form (cf. Appendix 0), and using different 
functional forms for the factor multiplying x a (1 — x) b . We have not seen any dependence of 
the uncertainty estimates on these changes. Although more radical ways of parametrizing 
the nonperturbative PDFs might affect the result more, there is no known example of such 
a parametrization, which at the same time still provides an equally good fit to the full data 
set. 

5 Further Examples 

5.1 W ± production at the LHC 

A study similar to the last section has been carried out for inclusive W ± production at the 
LHC. Figure H shows Xgiobai versus &w for the process pp — > W ± X at y/s = 14TeV, summed 

W production at the LHC 
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Figure 8: Minimum Xgiobai versus &w in n b, for inclusive W ± production at the LHC. (Cf. 
Fig. |3|.) The prediction is 189.7 nb. The points are the results of LM calculations. The curve 
is a polynomial fit to the points. 

over the two final states. The curve is a smooth interpolation of a series of PDF sets {S a } 
generated by the Lagrange Multiplier method. The best estimate value of the LHC cross 
section is aw = 189.7 nb. 

Comparing Figs. |3| and |8|, one immediately notices that the uncertainty of an/ (LHC) 
is greater than that of criy(Tevatron) for the same AXg lobal . This indicates that because W 
production in pp collisions at the LHC and pp collisions at the Tevatron involve different 
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mixtures of parton subprocesses as well as different kinematic ranges, the constraints imposed 
by current experiments included in the global analysis are also different for the two cases. 
Referring to the map of the d- dimensional PDF parameter space on the left side of Fig. [I], 
we are generating sample PDFs along different directions Lx in the two cases. Therefore it 
is not surprising that the rate of variation of Xgiobai i s also different. 

To demonstrate this point, and to quantify the uncertainty on the LHC prediction, we 
have carried out the same error analysis as in Sec.|], i.e., comparing the alternative PDFs 
to the individual experiments. Figure |^ gives the final overview of the 90% CL ranges of 
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Figure 9: Same as Fig. [?|, except for the LHC case. 



aw obtained from these comparisons, analogous to Fig. [7] for the Tevatron cross section. 
There are some differences compared to the Tevatron case. The LHC prediction is more 
tightly constrained by experiments that are sensitive to PDFs at small x. This makes sense, 
because W production at the LHC is not dominated by valence qq interactions. We note in 
particular that the two inclusive jet production experiments place significant constraints on 
aw at the LHC. 

We can combine the individual error bars in Fig. [5] according to the algorithm proposed 
in Sec. f| to produce a global uncertainty measure for aw at the LHC. The lower bound on 
aw0->HC) obtained by the intersections of the individual ranges is aw = 175.3 nb; the upper 
bound is aw = 204.6 nb. These bounds, shown as the dashed lines in Fig. [], correspond to 
±8% deviations from the prediction (189.7 nb). The global uncertainty on o"i4/(LHC) is thus 
significantly larger than that on aw (Tevatron). Reinforcing this conclusion is the fact that 
the scatter of the points in Fig. ^| is larger than in Fig. [7]. 

We can again inspect the increase in Xgiobai fr° m the quoted range of alternative fits. A 
±8% deviation from the minimum, symmetrized for simplicity, corresponds to the increase 

2 00 ^ 1 t t < ■ 1 mil il 1 at 1 . ' i m ml .it \~ i \ ill,-, i i i t , \ . i < • / i 1 11 \ - — 

of a w at the Tevatron. 



This number is similar to the increase in Xgiobai f° r our estimated uncertainty 



In the companion paper fllPf , we make some process- independent estimates of AXg lobal 
based on completely different considerations. Those arguments also yield the same order- 
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of-magnitude estimates of AXg lobal (in the range from 100 to 200) for acceptable PDFs 
around the global minimum. Since the effective Xgiobai> as a measure of goodness-of-fit, does 
not have a normal statistical implication, points on a constant Xgiobai surface in the PDF 
parameter space do not necessarily correspond to a constant likelihood. Some variation with 
the direction in the multi-dimensional space is to be expected. 

5.2 Uncertainties on Z° production 

We conclude this section by presenting results from applying the Lagrange Multiplier method 
to Z° production at the Tevatron and the LHC. 

Figure p!U| a shows the minimum Xgiobai as a function of oz at the Tevatron. The global 
prediction is oz = 6.55 nb. The experimental measurement by the DO collaboration is 
a z B = 0.221 ±0.003±0.011nb; the result from CDF (all data from Run I) is a z B = 0.250± 
0.004±0.010nb. (Here B is the branching ratio for Z° -> ee, which is (3.367±0.005) x 10~ 2 .) 
The comparison of the prediction to Tevatron data is discussed below. Analyzing the local 
Xn m the manner of Sec. || in order to quantify the uncertainty of the prediction, we find that 
the uncertainty of cr^ (Tevatron) due to PDFs is ±3% of the prediction. The corresponding 
increase in Xgiobai> symmetrized for simplicity, is approximately 130. 
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Figure 10: Minimum Xgiobai versus <7z for: (a) the Tevatron; and (b) the LHC. (Cf. Fig. |3|.) 
The predictions are 6.55 nb at the Tevatron and 58.0 nb at the LHC. The points are results 
of LM calculations. The curves are polynomial fits to the points. 

For the LHC process pp — > Z°X, Fig. |lO| b shows the minimum Xgiobai as a function of 
az- The dependence of Xgiobai on ^wy^LHC) exhibits a behavior departing from quadratic 
over the full range of aw/z under study. This is evidence that the Lagrange multiplier 
method can go beyond the traditional error matrix approach (which depends on the quadratic 
approximation) in exploring the neighborhood of the minimum. 

The global prediction is oz = 58.0 nb. Analyzing the local Xn as iu the other cases, 
we find that the uncertainty of cxz(LHC) due to PDFs is approximately ±10%. As in the 
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case of W ± production, the PDF uncertainty for Z° production at the LHC is significantly 
larger than that at the Tevatron. Measurement of and Z° production at the LHC will 
therefore provide significant information on PDFs. 



5.3 Comparison with existing data 
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Figure 11: Experimental measurements of ay/Bw and azBz at the Tevatron, compared to 
the PDF uncertainty band of the theoretical prediction. On the data points, the short error 
bar is the statistical error, and the long error bar is the combined statistical and systematic 
error. 

For W and Z production at the Tevatron, we can compare our calculated cross sections 
aw and az, with their ranges of uncertainty ±4% and ±3% respectively, to the measurements 



of CDF and DO from Run I. |26| The comparison is shown in Fig. [U]. The two experiments 
do not measure aw and az per se, but rather awBw and azBz where Bw is the branching 
ratio for W~ — > ev and Bz is the branching ratio for Z° — > ee. We have used the values 
Bw = 0.106 and B z = 0.0337 for the calculations [p7fl . The bands in Fig. [U] show the ranges 
of awBw and azB z from our PDF uncertainty study (but no uncertainty included from B w 
and Bz). The two measurements of awBw are consistent with the uncertainty range. The 
two measurements of azBz are not. 

It should be noted that CDF and DO use different normalizations for their luminosity 
determinations. The CDF collaboration bases its luminosity purely on its own measurement 
of the inelastic pp cross section [|28|, |29|, while DO uses the world average for this cross 
section. Thus current luminosities quoted by CDF are 6.2% lower than those quoted by DO. 
Consequently, all CDF cross section measurements are ab initio 6.2% higher than those of 
DO. If the CDF/DO measurements of awBw and azB z are rescaled by 6.2% with respect 
to each other, they are in excellent agreement. 

Because of the uncertainty in the inelastic pp cross section, it has been proposed to 
normalize future Tevatron (and LHC) physics cross sections to the measured W cross section 
(or rather awBw)- This makes the determination of the uncertainty of aw due to PDFs 
even more important. 
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5.4 Comparison with previous uncertainty estimates 

It is interesting to contrast our results to existing estimates of the uncertainties of aw and 
a z at the Tevatron and LHC colliders based on the traditional method of comparing results 
obtained from somewhat ad hoc PDFs. Some of these previous comparisons for aw (Tevatron) 
between historical PDFs as well as various trial up/down sets obtained by the CTEQ and 
MRST groups were shown in Fig. |2|. We will briefly comment on the results of |13[ in the 
context of this paper. 

Reference [[13] constructs an extended set of MRST PDFs, of which the most important 
for a w and az are the standard set MRST99 and three pairs of up/down sets designated 

{as has I}, {9 T,5 I}, {<? T, <? 1} 

in which some aspect of the parton distributions is either raised (f) or lowered (j) by an 
amount that represents an educated guess of a "standard deviation" . The predictions of aw 
and az are then compared for these alternative PDF sets to get an idea of the uncertainty 
due to PDFs. 

In the case of the Tevatron processes, the deviations of aw or az from the value for 
MRST99 for sets {a s T,"s |}, {g ],g !}, {q T,9 1} were found to be ±2%, ±1%, ±3% 
respectively. From these results, the authors of |13| concluded that the uncertainties of 
a w and az at the Tevatron are no more than about ±4%, and mainly attributable to the 
normalization uncertainty in the input u and d distributions. This conclusion appears to be 
quite consistent with the results of the previous sections based on exploring the variation 
of the cross section over the entire PDF parameter space. (This range of uncertainty also 
happens to coincide with what one would get by comparing historical PDF sets, as shown 
in the right plot of Fig. |^.) 

For the LHC, the MRST study found that the uncertainty of aw and az at the LHC 
is only slightly larger than at the Tevatron; the uncertainty was estimated to be ±5%. The 
largest observed variations came from the sets as f and as j, differing from the standard 
prediction by ±4-5%. This estimate is considerably lower than the ±8-10% result obtained 
by the detailed analysis of the previous sections. We have verified that the PDF sets that give 
±8-10% deviations of cxy^ (LHC) and az (LHC) from the standard prediction (represented by 
the points at the outer edges of the corresponding plots in Figs. || and [TO]) provide equally 
good or better fits to the global data sets compared to the fits of {as T> as Thus, 
it is clear that the Lagrange Multiplier method can generate optimal PDFs, i.e., having 

k The values of X^iobai f° r { a S 1\ a S 1} are {1531, 1356} compared to ~ 1400 for the outermost LHC 
sets shown in Figs. |j and [l^ and presented in the next section. 

'The global fits used at first in our exploration of PDF uncertainty were conducted with fixed a s . To 
make sure that this restriction does not result in an underestimate of the uncertainties of aw and az , we 
have examined the effect of freeing a s in the analysis (but imposing the known constraints from the world 
average of a s ). The results on the size of the uncertainties are not changed noticeably. This is because the 
full variations in the PDFs (particularly the gluon) allowed in the Lagrange approach can absorb the added 
degree of freedom. 
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the largest excursion of the variable X of interest, which are difficult to discover by ad hoc 
trial- and-error methods used in the past. 



6 PDF Sets for exploring W and Z Physics 

The parton distribution sets used in the above calculations are useful for exploring some 
aspects of W and Z physics at the Tevatron and LHC, since they provide much more re- 
liable estimates on the PDF uncertainties than existing ones in the literature, which are 
not designed to probe the full range of possibilities in the parton parameter space. With 
this in mind, we present in this section some representative PDF sets for applications to 
the rate of W and Z production at the Tevatron and LHC. These PDFs are relevant to the 
total cross sections aw and o~z, and each corresponds to a particular direction (Lx) in the 
PDF parameter space (see Fig. [I]). Therefore they are not suitable for estimating the PDF 
uncertainties of other observables that are sensitive to other aspects of the PDFs. Other 
PDF sets can be obtained, using the method introduced in this paper, to probe the range 
of other variables, such as rapidity (y) or transverse momentum (p?) distributions (hence 
relevant to the measurement of W mass). These will be investigated in subsequent work. 
Also, the companion paper |T(J supplies information from the Hessian method that can be 



used to construct the optimal PDFs for any observable X. 

The PDF set that yields the "best estimate" for all of the physical cross sections covered 
in this paper is our standard set Sq. The parametrization of the initial distribution is given 
in Appendix |Cj. In the following, we present two sets of PDFs that bound the likely range 
for each of the cross sections. 

To exemplify the PDFs that characterize the range of uncertainty of W production 
at the Tevatron, we use two representative sets, labeled S WTeV , which correspond to aw = 
cw{So)±Aaw (with Aaw/^w ~ 0.04) respectively. These two sets are extreme fits obtained 
by the Lagrange Multiplier method. The parameters {a} for these sets are given in Appendix 
0. We now compare some of the parton distributions from the three sets (Syy x e v ^o, iw)> 
to examine the ranges of variation of the PDFs themselves. 

Figure |12| a shows u(x,Q), d(x,Q) and g(x,Q) for S WTeV , compared to the standard 
set S , at Q = 80 GeV. The function xAf(x) is plotted for each parton flavor, where Af is 
f — fo- The gluon function has been divided by 10 to fit on the same graph. The solid curves 
(Au+, Ad + , Ag + ) correspond to S WTeV , and the dashed curves to S WTeV . For S WTeV , 
requiring aw to be larger than a w makes the u and d distributions larger than for the 
standard fit (w and d ) so Au + and Ad + are positive. Then the gluon distribution must 
be smaller than the standard because of the momentum sum rule. In the case of S W TeV , 
the reverse is true, resulting in almost a mirror behavior. At the Tevatron, a typical x for 
the parton-level process q{q 2 — > W is Mw/yfs = 0.04. The differences Am and Ad are 
significant in the range 0.01 < x < 0.04. The magnitude of Af(x) in this range is a few 
percent of the standard fo(x), which makes sense since Aaw is a 4% shift of aw for these 
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Figure 12: Comparison of PDF sets S w to the standard set So for the Tevatron (upper) and 
LHC (lower) cases respectively. Au± is the difference between u(x) with aw = <?wo ± Aaw 
and u{x) with aw = fivo! Q = 80 GeV. The solid curves are xAf + (x) and the dashed curves 
are xAf~(x). The abcissa is \og 10 x. 



PDF sets. 

We can carry out the same comparison for W production at the LHC. The PDFs 
that bound the range of uncertainty are designated as S WLUC , which correspond to aw = 
vwiSo) =t Aaw with Aaw/o'w ~ 0.08. The PDF parameters are given in Appendix [Cj. 
Figure |l2"l b shows parton distributions from S WLHC . (Again, the gluon has been divided by 
10.) In the LHC case, the typical x for the process qiq2 — > W ± is Mw/y/s = 0.006. The 
region where Am and Ad are significant is seen accordingly lower in x than for the Tevatron 
case. 

In Appendix |C| we also present PDF sets S% TeV and S% LHC that characterize the range 
of uncertainties of Z production at the Tevatron and LHC. These correspond to the outlying 
points on Figs. |l0|a,b. They are similar to S W TeY and 5^ LHC , with small differences in the 
flavor dependence. 
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7 Summary 



We have developed the Lagrange Multiplier method to calculate uncertainties of physical 
observables due to PDFs, and we have used the method to determine the uncertainty of 
the total cross sections for W production and Z° production at the Tevatron and LHC. 
The method is more reliable than past estimates because: (i) it explores all the possibilities 
in the parameter space of the input PDFs, independent of other assumptions; and (ii) it 
produces the maximum allowed range for the specified physical variables. This is in contrast 
to previous attempts which relied on varying certain features of the parton distributions 
chosen in some ad hoc way. 

From this analysis, we find that the uncertainty of the prediction for aw or oz at 
the Tevatron with current experimental constraints is approximately ±3-4%; and at the 
LHC the uncertainties are approximately ±8-10%. These numbers do not include other 
uncertainties associated with theoretical approximations, nuclear corrections, and other un- 
expected sources. We have explored to some extent the possible effects due to the choice of 
parametrization of the nonperturbative input PDFs, and found them to be small. The cur- 
rent work should be considered exploratory in nature, as a first application of this improved 
approach to error estimates. A more comprehensive study, based on soon to be improved 
data sets, and including other sources of uncertainties, will produce better overall estimates 
of the physical predictions. 

This study should be regarded as the precursor for many interesting applications to 
come, on physical processes of interest to the precision study of the Standard Model, and on 
predictions for New Physics at future colliders. Some examples are rapidity distributions of 
W and Z° production, which contain a wealth of information on parton structure of the 
nucleon; the W mass measurement; top and Higgs cross sections, etc. 

There are other approaches to error estimates in global QCD analysis |], 0, 0, ||. In 
general, if greater emphasis is placed on the "rigor" of the statistical method, then the range 
of experiments that can be included in the analysis is narrower. We have chosen to emphasize 
the inclusion of the full range of experimental constraints, and adapt the statistical analysis 
to deal with the practical problems that one faces in such a system. Within our general 
framework, there is an alternative, complementary approach based on the conventional error 
matrix method fl. We explore this latter method, as applied to global QCD analysis of 



PDFs, in a companion paper fllQ[ . We mention briefly the contrasting features and relative 
merits of the two approaches here. 

The Lagrange Multiplier method focuses on a given physical observable X (or a set 
of observables {X k }) and determines the uncertainty AX allowed by the global data set 
within a specified tolerance for the global fit. The error matrix approach, using the Hessian 
matrix, focuses instead on the uncertainties of the PDFs as represented by the parameters 
{cii] i = 1, . . . , d}. It is in principle universal because, once determined, these errors can be 
propagated to any physical variable X. However, the results are reliable only if the function 
Xgiobai( a ) an d the observable X(a) can be approximated by lowest order expansions in the 
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parameters {a}, and if the numerical computation of the derivatives (the Hessian matrix) is 
under control. The latter problem is surprisingly difficult for global QCD analysis, because 
the eigenvalues of the error matrix vary by many orders of magnitude. This problem has 
been solved f|, and the error matrix results are consistent with the constrained fitting results 
TCfl . Thus, at present, both methods appear to be applicable to the study of uncertainties 



in global QCD analysis. 

In Figures || and |10| b there is a significant cubic term in the dependence of Xgiobai on 
<7vi/(LHC) and <Tz(LHC), respectively. To calculate Xgiobai versus X accurately in such cases, 
the Lagrange Multiplier method is necessary. Traditional linear error analysis based on the 
Hessian matrix can only produce a quadratic approximation to the dependence. 

When both methods are applicable, the Hessian method is more flexible and easier to 
apply computationally. But generally the Lagrange method is more robust and reliable. As 
we expand the investigation to other physical processes of interest, we will continue to test 
the efficacy of both methods and cross check the results. 
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A The effect of correlated errors on A% 2 

The global fitting function Xgiobai defined in (|2|) resembles the standard statistical variable x 2 , 
so it is tempting to try to apply theorems of Gaussian statistics to analyze the significance 
of the fit between theory and experiment. However, the familiar theorems do not apply, 
because of correlations between measurement errors. The purpose of this Appendix is to 
explore this issue. The effect of correlated errors is potentially a source of confusion. 

For simplicity we describe the simplest case: measurement of a single observable. The 
arguments can be extended to cases where multiple quantities are measured, such as the 
determination of parton distribution functions. 

Consider an observable m that is measured N times. We shall refer to N measure- 
ments of m as one "experiment". Let the true value of m be m . The measurements 
are m%, m 2 , m 3 , . . . , m N . The deviations from the true value are «i, a 2 , a 3 , . . . , a^, where 
«i = rrti — itlq. In general the measurement errors are correlated, so in the Gaussian approx- 
imation the probability distribution of the fluctuations is 

f 1 N } 

dP = Afexp < — ctiCijaj > d N a. (11) 

I t,i=i J 

Here is a real symmetric matrix, and M = VDet C / (2tt) n ^ 2 ensures the normalization 
condition J dP = 1. 
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We will need the variance matrix (a^a^), where the notation (Q) means the average of 
Q in the probability distribution QllD . For this Gaussian distribution, 



<o i a i > = (C7-%. (12) 

The mean square fluctuation Ei of the 2 th measurement is 

E i =(a& = (C- 1 ) u . (13) 

To find the best estimate of the value of m from these N measurements, ignoring the corre- 
lations in the measurement errors, we define a chi-squared function Xu{ m ) by 



.2, 



N ( ,2 

[rrii — m) 



E oii — in i . . 

E , • ( 14 ) 

i=l * 

The value of m that minimizes Xu( m )i can ^ i s then the best estimate of m based on this 
information. The function Xu( m ) is analogous to the fitting function Xgiobai i n ^ n e CTEQ 
program, in the sense that it does not include information about the correlations between 
errors. The minimum of Xu( m ) occurs at a weighted average of the measurements, 

If all the EiS are equal then m is just the average of the measurements. 

Now, what are the fluctuations of the mean ml That is, if the "experiment" consisting 
of iV measurements could be replicated many times, what would be the distribution of m's 
obtained in those many trials? It turns out that m has a Gaussian distribution 

dP 1 

* = 7i^ expHra - mo)2/(2S2)1 - (16) 

The standard deviation S of m is the RMS fluctuation; that is, 

in 3 



where 

D = y- 

^ E t 



The question we wish to answer is this: How much does Xu( m ) increase, when m 
moves away from the minimum (atfn) by the amount ±E that corresponds to one standard 
deviation of the mean ? The answer to this question is 

Axt = T?D. (19) 
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This result follows easily from the definition (fH]), because 

Xl(m + E) - xl(m) = -2E £ + E 2 £ i-, (20) 

i i 

and the term linear in E is by the definition of m. So far the discussion has been quite 
general. We will now examine some illustrative special cases. 

Example 1: Suppose the measurement errors are uncorrelated; that is, 

dj = Sij/Ei. (21) 



Then the standard deviation of the mean m is E = 1/yD. Thus for the uncorrelated case, 
the increase of x\ corresponding to one standard deviation of the mean is Ax 2 = 1. This is 
the "normal" statistical result: The la range corresponds to an increase of x 2 by 1. 

An even more special case is when the errors are uncorrelated and constant: Ei = a 2 
independent of i, where a is the standard deviation of single measurements. The correlation 
matrix is Cij = 5ij/a 2 . In this case D is N/a 2 , and the standard deviation of the mean is 
E = a/y/N. 

The criterion Ax 2 = 1 for one standard deviation of a measured quantity is a standard 
result, often used in the analysis of precision data. But if x 2 is defined ignoring the correla- 
tions between measurement errors, then the criterion A% 2 = 1 is only valid for uncorrelated 
errors. We will next consider two examples with correlated errors, to show that A% 2 is not 
1 for such cases. 

Example 2: Suppose measurements 1 and 2 are correlated, 3 and 4 are correlated, 5 
and 6 are correlated, etc. Then the correlation matrix is 

1 jo 2 for % = j 

C i:j = { c/a 2 for ij = 12 or 21, 34 or 43, etc (22) 
otherwise 

where — 1 < c < 1 since the determinant of C must be positive. The inverse matrix C~~ l can 
be constructed using the fact that C is block diagonal, consisting of N/2 2x2 blocks. Then 
it can be shown that 

and Axl = 1 - c. (23) 



VN^/TTc 

The increase of x\ f° r one standard deviation of the mean ranges from to 2, depending on 
c. The criterion A% 2 = 1 does not apply to this example with correlated errors. A standard 
increase of x\ ma y be smaller or larger than 1. 

Example 3: For an even more striking example, suppose the N measurements that 
constitute a single "experiment" are, for i — 1, 2, 3, ... , N, 

rrii = m + i/i + P (24) 
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where the yi are randomly distributed with standard deviation cr, and the measurements 
are systematically off by the amount ft. Suppose that ft has a Gaussian distribution with 
standard deviation s for replications of the "experiment" . In this example, 

c « = M^-jto)' <25) 

(C-% = v'Stj + s 2 . (26) 
The variance of the individual measurements (rrii) is 

(m 2 ) - (m) 2 = cr 2 + s 2 . (27) 
Therefore our uncorrelated chi-squared variable X«( TO )> defined ignoring the correlations, is 

N . ,2 

Xl(m) = J2 2 ■ 2 • ( 28 ) 
z — ' cr z + s 

i=i 

The minimum of x«( m ) occurs at m, which is just the average of the individual measure- 
ments. The variance of fa, averaged over many replications of the "experiment" , is 

2 

s 2 = (m 2 ) - (m) 2 = s 2 + (29) 

The increase of x\ as m moves from m to fn ± S, ie., by one standard deviation of the 



a 2 + Ns 2 
a 2 + s 2 



= XZ(rn + E) - = 2 , 2 ■ (30) 



In the limit s/a 1, the error correlations in this model become negligible and Ax 2 reduces 
to the conventional value of 1. But in the limit s/a ^ 1 where the error correlations are 
dominant, A% 2 approaches N. 

Thus for Example 3 — a systematic error with 100% correlation between measurements — 
the increase of xt f° r a standard deviation of m is much larger than 1. If s and a are 
comparable, then A^ 2 is of order N. 

If the correlation matrix CV, is known accurately, then the correlation information can 
be incorporated into the definition of the x 2 function, in the manner of Appendix |B|. For 
the full list of experiments in the global analysis of parton distribution functions, however, 
the correlations of systematic errors have not been published, so the fitting function Xgiobai 
has only uncorrelated systematic errors. 

We described above the measurement of a single quantity. The determination of par- 
ton distribution functions seeks to measure many quantities, i.e., the 16 parameters {a}. 
The above arguments can be extended to measurements of multiple quantities. If the mea- 
surement errors are uncorrelated, then the increase of x\ by 1 from the minimum defines a 
hyperellipse in parameter space — the error ellipse — corresponding to one standard deviation 
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of linear combinations of the parameters. However, if the errors are correlated then Axi — 1 
is not the correct criterion for a standard deviation. 

The Lagrange Multiplier method finds the best fit to the data, subject to a constrained 
value of some quantity X. The prediction of X is at the absolute minimum of xt- Again, 
if the errors are uncorrelated then one standard deviation of the constrained quantity corre- 
sponds to an increase of x 2 by 1 from the absolute minimum. But if the errors are correlated 
then Axi = 1 is not the correct criterion for one standard deviation of X. 

One reason for describing this familiar, even elementary, statistics, is to avoid certain 
misconceptions. Our standard PDF set So is a parametrized fit to 1295 data points with 16 
fitting parameters. The minimum value of Xgiobai * s approximately 1200. Naively, it seems 
that an increase of Xgiobai by merely 1, say from 1200 to 1201, could not possibly represent 
a standard deviation of the fit. Naively one might suppose that a standard deviation would 
have Ax 2 ~ vT295 rather than I. However, this is a misconception. If the errors are 
uncorrelated (or if the correlations are incorporated into x 2 ) then indeed Ax 2 = 1 would 
represent a standard deviation. But this theorem is irrelevant to our problem, because the 
large correlations of systematic errors are not taken into account in Xgiobai- 



B x 2 function including correlated systematic errors 

The purpose of this appendix is to derive the appropriate definition of x 2 f° r data with 
correlated systematic errors. The defining condition is that x 2 should obey a chi-squared 
distribution. 

Let {mi} be a set of measurements, where i = 1, 2, 3, . . . , N. Let tj be the true, i.e., 
theoretical value of the i th measured quantity. Several kinds of measurement errors will 
contribute to the difference between rrii and U. The uncorrelated error of measurement i 
is denoted by Oi. There are also correlated errors, K in number, denoted /3u, foi, ■ ■ ■ , Pxi- 
Thus the i th measurement can be written as 

K 

mi = U+ errors = U + + ^ (3^ (31) 

where and r'- are independently fluctuating variables. We assume that each of these 
fluctuations has a Gaussian distribution with width 1, 

P(r) = e -^. (32) 



2,71 

Note that r' is independent of i; that is, the errors /3ji, /3j2, ■ ■ ■ ,PjN are 100% correlated for 
all iV data points. 

The probability distribution of the measurements is 

N K 



dP = Y^pir^dnY^pir'^dr'j 

J i=i j=i 
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N 



K 



X 



rrii 



(33) 



Now we will evaluate the integrals over and r'- in two steps. First evaluate the rj integrals 
using the delta functions, 



dP= / Yldr'f.e-^ 2 d N m 

J .7=1 



(34) 



where Ci is a normalization constant and 

X ' x ' 



:? = E 



i=l 



(35) 



Note that x? is a function of rj,...,r^. These variables {r'j} could be used as fitting 
parameters to account for the systematic errors: Minimizing xl with respect to r'- would 
provide the best model to correct for the systematic error of type j. Because xl is only a 
quadratic polynomial in the r'- variables, the minimization can be done analytically. 

To continue evaluating ( |33| ) we now do the integration over {/•'•}. Write m the form 



N , ,2 K K 

t=i 4 i=i jj'=i 



(36) 



where -Bj is a vector with K components 



N 



Bj = ^Pjiirrii - U)/a 



(37) 



and Ajji is a K x K matrix 



N 



i=l 



Then the integration over d r' is an exercise in Gaussian integration, with the result 



dP = C exp 



where C is a normalization constant and 



2 X 



d N m 



(38) 



(39) 



N , , \ 9 K K 

Tfli — Is) 



X 



E 

i=i 



EE" •>..,«•• 

3=1 j'=l 



(40) 



This equation is the appropriate definition of x 2 for data with correlated systematic errors. 
The correlated errors are defined by the coefficients /3ji in ([H]), which determine the vector 
Bj and matrix Ajj>. An interesting relation is that the x 2 quantity in ( f40|) is the minimum 
of xl with respect to the parameters r[, . . . ,r' K . 
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Another expression for % 2 , which may be derived from by Gaussian integration, 



is 



N N 



X' 



(41) 



i=l i'=l 



where Va is the variance matrix 



K 



3=1 



(42) 



It can be shown that the inverse of the variance matrix is 

3,j'=l 



(43) 



Therefore (HO) and (HI) are equivalent. However, there is a real computational advantage in 



the use of (IT) because it does not require the numerical inversion of the N x N variance 
matrix. 

To check that (|4~0D makes sense we can consider a special case. Suppose the number K 
of systematic errors is N, and each systematic error contributes to just one measurement. 
Then the matrix of systematic errors has the form 







(44) 



This situation is equivalent to an additional set of uncorrelated errors {&«}. The vector Bj is 
then 



and the matrix Ajf is 



B, 



Ajj' — Ojj> 



O": 



1 + 



Substituting these results into fl40|) we find 



2 _ (m, - tj) 



(45) 



(46) 



(47) 



which makes sense: the uncorrelated errors just combine in quadrature. 

The statistical quantity y 2 has a chi-squared distribution with iV degrees of freedom. 
Thus this variable may be used to set confidence levels of the theory for the given data. But to 
use this variable, the measurement errors o"j and j3ji, for i = 1, 2, . . . , iV and j — 1, 2, . . . , K, 
must be known from the experiment. A chi-squared distribution with many degrees of 
freedom is a very narrow distribution, sharply peaked at y 2 = N. Therefore small inaccuracies 
in the values of the aVs and fy's may translate into a large error on the confidence levels 
computed from the chi-squared distribution. 
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It is equation (pEOD that we use in Section |] to compare the constrained fits produced by 
the Lagrange multiplier method to data from the HI and BCDMS experiments. Correlated 
systematic errors are also used to calculate \ 2 f°r the CDF and DO jet experiments. 



C Parton Distribution Sets 

We give here the PDFs described in Section 6. So is the standard set, defined by the 
absolute minimum of Xgiobai- $w,Tev are ^s to the global data sets with extreme values of 
<Tvi/(Tevatron), i.e., the outermost points on Fig. [| generated by the Lagrange multiplier 
method. Sg Tev , S^ LHC , and S% LHC are analogous for Z production and W and Z production 
at the LHC. 

The functional form of the initial parton distributions and the definitions of the PDF 
parameters at the low-energy scale Qq — 1 GeV are 



f{x,Ql) 



A x Al (l-x) A2 (l + A 3 x M ) 



for u v , d v , g,u + d, s(= s); and for the ratio 

djx, go) 
u(x,Ql) 

The tables of coefficients follow. 



X) 



+ (1 + A 3 x)(l 



x 



.A 4 



S : 





A 


A x 


A 2 


A 3 


A, 




0.5959 


0.4942 


4.2785 


8.4187 


0.7867 


Uv 


0.9783 


0.4942 


3.3705 


10.0012 


0.8571 


9 


3.3862 


0.2610 


3.4795 


-0.9653 


1. 


d/u 


3.051E4 


5.4143 


15. 


9.8535 


4.3558 


u + d 


0.5089 


0.0877 


7.7482 


3.3890 


1. 


s 


0.1018 


0.0877 


7.7482 


3.3890 


1. 
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A 

Aq 


A 

M 








d v 


0.2891 


0.5141 


3.8555 


in n c on 

10.9580 


0.4128 




0.2184 


0.2958 


4.6267 


35.7229 


1.0958 


u v 


1 A1 ,1 O 

1.0142 


0.5141 


3.3614 


n onn r 

9.2995 


n o n r o 

0.8053 




r\ c\r\ 

0.2979 


0.2958 


O Of* ^7(^\ 

3.3279 


32.8453 


0.9427 


9 


4.6245 


r\ A or A 

0.4354 


3.4795 


—0.9728 


1. 




1.8080 


0.0458 


3.4795 


—0.0519 


1. 


d/u 


0.9UO-C/4 


5.66 ^ 3 


15. 


9.o535 


4. f45o 




2.041E4 


5.1506 


15. 


9.8535 


4.8320 


u + d 


0.4615 


0.0108 


6.6145 


0.92784 


1. 




1.2515 


0.3338 


7.5216 


-0.0570 


1. 


s 


0.0923 


0.0108 


6.6145 


0.9278 


1. 




0.2503 


0.3338 


7.5216 


-0.0570 


1. 





A 


A 1 


A 2 


^3 


A, 


d v 


0.6061 


0.5502 


4.0017 


5.8346 


0.5343 




0.3427 


0.3728 


4.5166 


19.8510 


0.9966 


u v 


1.2159 


0.5502 


3.3347 


7.3386 


0.7711 




0.5247 


0.3728 


3.3905 


20.1006 


0.9556 


9 


4.4962 


0.4321 


3.4795 


-0.9023 


1. 




2.3113 


0.1032 


3.4795 


-0.6349 


1. 


d/u 


4.321E4 


5.4724 


15. 


9.8535 


4.6298 




2.818E4 


5.4540 


15. 


9.8535 


4.4376 


u + d 


0.4609 


0.0103 


6.6671 


0.9822 


1. 




0.9900 


0.2926 


8.3205 


2.1648 


1. 


s 


0.0921 


0.0103 


6.6671 


0.9823 


1. 




0.1980 


0.2926 


8.3205 


2.1648 


1. 
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A 

A 


A 

M 


A 2 




A 4 


d v 


0.762b 


0.5008 


4.6393 


10.8532 


1.0595 




0.5671 


0.4771 


4.2615 


8.8355 


0.8130 


u v 


1.0608 


0.5008 


3.4023 


9.6622 


0.8968 




0.9142 


0.4771 


3.3761 


10.9138 


0.8809 


9 


2.2679 


0.0733 


3.4795 


—0.9860 


1. 




2.5021 


0.3981 


3.4795 


1.6229 


1. 


a/u 


O 1 *70 TP A 
2.1 (Qh/4: 


0.2o f 


Id. 


n o c o c 


a /loin 
4.481U 




4.531E4 


5.4979 


15. 


9.8535 


4.6585 


u + d 


1.1980 


0.2952 


6.9475 


-0.5442 


1. 




0.2759 


-0.0918 


8.2045 


6.3950 


1. 


s 


0.2396 


0.2952 


6.9475 


-0.5442 


1. 




0.0552 


-0.0918 


8.2045 


6.3950 


1. 





A 


A 1 


A 2 


^3 


A 4 


d v 


0.5659 


0.4616 


4.5297 


12.3685 


0.9836 




0.4585 


0.4496 


4.2122 


10.3850 


0.7760 


u v 


0.8344 


0.4616 


3.3847 


12.1129 


0.8872 




0.7640 


0.4496 


3.3566 


12.8253 


0.8701 


9 


2.3282 


0.0918 


3.4795 


-0.9837 


1. 




2.9475 


0.4219 


3.4795 


0.9447 


1. 


d/u 


2.421£4 


5.3032 


15. 


9.8535 


4.5341 




4.416E4 


5.4708 


15. 


9.8535 


4.7925 


u + d 


1.1130 


0.2698 


6.8490 


-0.5330 


1. 




0.2719 


-0.0899 


8.1492 


6.5300 


1. 


s 


0.2226 


0.2698 


6.8490 


-0.5330 


1. 




0.0544 


-0.0899 


8.1492 


6.5300 


1. 



With a program to solve the PDF evolution equations, the PDFs for an arbitrary momentum 
scale Q can be generated. 
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